Molecular Dynamics on Diffusive Time Scales from the Phase Field Crystal Equation 
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We extend the phase field crystal model to accommodate exact atomic configurations and vacan- 
cies by requiring the order parameter to be non-negative. The resulting theory dictates the number 
of atoms and describes the motion of each of them. By solving the dynamical equation of the model, 
which is a partial differential equation, we are essentially performing molecular dynamics simulations 
on diffusive time-scales. To illustrate this approach, we calculate the two-point correlation function 
of a liquid. 
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Molecular dynamics (MD) has long been a powerful 
tool to study statistical mechanical systems (for an in- 
troduction, see (e.g.) Ref. [1]). By postulating the inter- 
action between atoms and solving the resulting equations 
of motion, precise information about each atom is known. 
One of the drawbacks of MD, however, is that too much 
information is captured. For example, atomic motions 
in MD simulations are resolved on atomic time scales, 
whereas in many systems the relevant time scales are 
diffusive. This makes MD computationally demanding, 
if not completely inapplicable, in many cases of interest 
where long time scales are required. In this paper we 
pursue a novel approach to attain long time scales, start- 
ing not from individual particles but from a continuum 
description of matter known as the phase field crystal 
(PFC) modelg, 1,0,0,0. 

The starting point of the PFC model is that crystalline 
materials are governed by a free energy functional that 
penalizes departures from periodicity of the density in 
the same way that the Landau theory of phase transi- 
tions uses a functional that penalizes spatial gradients of 
the order parameter. The PFC model is formulated in 
terms of an order parameter representing the local den- 
sity, and is constructed so that the free energy functional 
is minimized by a periodic order parameter configura- 
tion. Despite its simplicity and minimal physical input, 
the PFC model can reproduce both qualitative and semi- 
quantitative (i.e. scaling) properties of multicrystalline 
solidification [4], dislocation dynamics[7], fracture, grain 
boundary energetics 0], elastic (phonon) interactions Q, 
grain coarsening 0, linear and nonlinear elasticity [4], and 
plasticity [To[. The PFC model has also been extended to 
binary systems (HI. and can be related to density func- 
tional theory [Bj]. Recent applications of the renormaliza- 



tion group technique 11 



[12|, [13[ and adaptive mesh re- 
finement have improved the computational efficiency of 
the model, with resultant computational times several 



Hi, [14]. 



orders of magnitude times faster than MD 

Although the PFC model represents microscopic con- 
figurations, it is not MD. The model describes the col- 



lective properties of the crystal, but it does not attempt 
to describe the motion of each individual atom. One 
can regard the peaks in the order parameter as repre- 
senting local density maxima, and thus be identified as 
PFC 'atoms'. However, although the order parameter, 
p(x, £), tends to form PFC 'atoms' in order to minimize 
the total energy of the system, their number is not con- 
served. This neglect of the actual atomic configuration, 
and the resulting absence of vacancies in the description, 
prevents us from using the model to describe faithfully 
microscopic phenomena that involve atomic hopping and 
vacancy diffusion. 

The goal of this paper is to modify the PFC model such 
that it describes not only the collective behavior, but also 
the motions of individual atoms. We will see that this can 
be done by constraining the value of the order parameter 
to be positive. By so doing, instead of being an abstract 
order parameter, p(x, t) becomes a physical density — the 
number of atoms in the model can be controlled by ad- 
justing a single parameter, po- The resulting theory is a 
MD simulation: we can specify the temperature, number 
of atoms and the interaction potential between atoms. 
As an illustration of this approach we simulate a simple 
liquid and reproduce the form of the standard two-point 
pair distribution function. 

Inclusion of Vacancies:- In real materials, vacancies are 
present when the local density is low, i.e., when there are 
not enough atoms to fill the space. In the PFC model, 
however, even if the value of the order parameter is small, 
which is analogous to the low density situation, a perfect 
periodic configuration can still be formed because there 
is no constraint, or energy penalty, for negative values of 
the order parameter. Therefore, as long as the system is 
in a periodic state, such as the 2-D triangular phase, any 
uniform configuration will evolve to a spatially periodic 
one in equilibrium. Thus, the notion of vacancies is not 
respected in this model. If a vacancy is created through a 
special initial condition, the free volume will simply dif- 
fuse throughout the crystal as the configuration readjusts 
its periodicity. 
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We can stabilize vacancies by imposing a constraint on 
the order parameter — we forbid the order parameter to 
be negative. In this case, if the local order parameter 
is not high enough, instead of forming a periodic state 
that extends to negative values, the system can form a 
periodic structure in some region, while leaving a very 
low, or zero, density in another. The number of atoms is 
then conserved and the zero density regions are identified 
with vacancies. 

We now identify the region of the phase diagram in 
which vacancies are present and stable, and we do this by 
calculating the energy of a state with vacancies, working 
for simplicity in two dimensions (2D). The PFC model is 
given by the free energy density [3|, [4J, 
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where r < is the undercooling parameter and po is the 
mean value of the order parameter. The dynamics follows 
the "Modified" PFC formulation!! 
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where F = J fd d x is the total free energy of the system, 
a and (3 are parameters that control the evolution, and r] 
is a Gaussian white noise satisfying the usual fluctuation- 
dissipation theorem. It is helpful to introduce the ansatz 
for the one-mode approximation to the triangular state 
in two-dimensions, 



p(x) = Aj2(( 



-ik. 



Po, 



(3) 



where £1,2,3 = (V3/2)y± (l/2)x are the basis wavevec- 
tors of the triangular phase. Substituting this ansatz into 
Eq. (pp), and averaging over the whole system gives the 
free energy density as a function of the constant ampli- 
tude, A: 



fo( P o, A) = —A i - 12A 3 p + ^(2 + 2r + p 2 ) 
+ 3A 2 (r + 3p 2 ). 



(4) 



Minimizing fo(po, A) with respect to A gives two roots 
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A±(po) = ^ i^Spo ± y/-15r - 36p§J , (5) 

where the solutions that minimize the energy are A = A+ 
for po > 0, and A = A- for po < 0. The roots are real 
for po < v / -5r/12 (recall that r < 0). 

Now, let us consider the effect of the constraint that 
the density be positive. Examining Eq. (|3]), we see that 
the summation is bounded by ±6, so requiring p(x, t) > 
is equivalent to requiring \A\ < po/6. However, Eq. ([5]) 



shows that |A + (po)| > Po/6 for all values of r and po, 
so the ground state A = A+ is forbidden by the con- 
straint. The ground state must be given by some other 
configuration. 

There are at least two possible configurations for the 
ground state. First, the ground state can still be perfectly 
periodic with an amplitude A ^ A+ satisfying \A\ < 
Po/6. Second, the ground state can partition itself into 
two domains — a perfectly periodic domain with average 
density pi and amplitude A\ satisfying \A\\ < pi/ 6, and 
a domain with p(x) = 0. The second domain corresponds 
to vacancies. To see which is realized in practice, we have 
to calculate the energy of these two states, and recognize 
that the ground state is the one with lower total energy. 

Let us first calculate the free energy density of a per- 
fectly triangular state. Since A = A + is forbidden, we 
are left with three options for A: A = A_, which is the 
other local minimum of the free energy, and A = ±po/6. 
The latter two are the boundary values satisfying the 
condition \A\ < po/6. By examining Eq. (|4]), one can 
see that /o(po>Po/6) < /o(po, -po/6), so we can ignore 
the A = — po/6 solution. The free energy density of the 
periodic state is then 



(6) 



fper(po) = fo (PO, ^) 

if |t4_(po)| > Po/6, and otherwise, 

f P er(po) = Min (/o (po, A-(po)) Jo (po, y)) (7) 

where Min(a, b) denotes the minimum of a and b. Sub- 
stituting A_(po) and po/6 into Eq. (j4j) gives the explicit 
expressions 



f (po,A-(p )) = -^Po 
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-15r-36p§ (8) 



/o ( p0 ' f) = 2^ [133 ^ + (M4 + 168r) ^ ] * (9) 

Now, let us compare the energy of these two possible 
ground states. If the system is perfectly periodic over the 
whole domain, whose area is designated Bo, then the free 
energy is given by 



fwhole(Po) — A)/per(Po), 



(10) 



If the whole system instead partitions itself into one do- 
main made up of a triangular phase having mean density 
Pi > Po, with the remaining domain having p = 0, the 
free energy is given by (for simplicity, surface energy be- 
tween the two phases is neglected in this calculation.) 

fv(po) = B 1 f per (pi) = Bofper(pl), (11) 

where B\ is the area of the triangular domain. The sec- 
ond equality is obtained by using the conservation of 
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mass PoBq = p\b\. The difference between these two 
free energies, Af = f v - f who i e , is 



A/ = BoPo 



fper(pi) _ fper(Po) \ 
Pi PO ) ' 



(12) 



It is important to note that pi is a parameter that we 
can choose to minimize the energy of the second possible 
state; the only constraint is that p\ > po because B\ < 
Bo. 
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FIG. 1: The function fo(po, A-(po))/ po is an increasing func- 
tion of po, plotted for various values of r. The x-axis is plotted 
on logarithmic scale in order to resolve the curves. 

For vacancies to exist, we require that Af < for 
some values of p\ > po. We note, however, that for 
the solution A = A_, /o(po> A-(po))/po is an increasing 
function of po (see Fig. (pQ)) and so Af is positive for 
this branch of the solution. In other words, no vacancy 
is present in this solution. Therefore, in order to have 
vacancies in the ground state, we require this branch of 
the solutions to be forbidden by the constraint; i.e., we 
require |A_(po)| > Po/6, which by Eq. |5|) is equivalent 
to requiring 



p ,pi < \/^12r/53. 



(13) 



On the other hand, it is easy to show from Eq. ([9]) that 
fo(po, Po/6)/po has a minimum (for r < —6/7) at 



^(-48- 560/133 



(14) 



Thus, if po < Pmin and r < —6/7, the system can mini- 
mize the total free energy by partitioning itself into two 
domains: a triangular phase made up of 'atoms', with av- 
erage density p — Pmin, and a region of vacancies where 
p = 0. Combining Eqs. ([T3j) and ([T4|) indicates that the 
minimum also satisfies the constraint |A_(po)| > Po/6 
so long as r > —636/343. We also note that the area 
of the triangular phase, B\ = Bo(po/ pi), is directly pro- 
portional to the mean density, po. So by adjusting po, 
we can control the number of atoms in the PFC model. 
This shows that the addition of the constraint, p(x) > 
for all x, does indeed promote the p(x) from an abstract 



order parameter to a physical density, which dictates the 
number of atoms in the system. 

To summarize, the various constraints define the region 



Pn 
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-48 - 56r 



133 



636 6 /irN 

and -_< r <-- (15) 



where the triangular phase and stable vacancies can coex- 
ist. The area fraction of the triangular phase is po/pmin, 
and the amplitude is A = pmin/6- The rest of the do- 
main has zero density and thus is composed of vacancies. 
These results are summarized in Fig. (|2]). 
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FIG. 2: Window in which stable vacancies can coexist with 
a triangular phase. A minimum in /o(po, po/6)/po exists at 
pmin — v (—48 — 56r)/133 for r < —6/7, and this min- 
imum is forbidden by the constraint |A_ (po) | < Po/6 for 
r > -636/343. 



Implementation:- In order to implement the positive den- 
sity constraint, we add a vacancy term, f vac (p), to the 
free energy density that penalizes negative values of 
p(x,t). As long as the repulsion from negative values 
is strong enough to avoid p < 0, the result should not 
depend on any particular choice of f vac (p)- Of the many 
possible choices for f vac (p), we use 



f vac (p) = H(\p\ n - p n ), 



(16) 



with n = 3 and H = 1500, because this turns out to be 
numerically convenient and stable. 

With the vacancy term, Eq. (fTS]). we can numerically 
verify the analytical calculation for the coexistence be- 
tween the periodic phase and vacancies. Fig. (|3|) provides 
results from simulations with r = —0.9 and different val- 
ues of po, showing clearly that the number of atoms in- 
creases with po- In addition, Fig. shows that the PFC 
atomic density (i.e., the number of atoms per unit area.) 
indeed increases linearly with po, as expected. The curve 
starts to saturate at around po = 0.15, as opposed to the 
prediction from Eq (fT5|), pmin ~ 0.134. This discrepancy 
is not surprising for several reasons: In the analytical cal- 
culation, we consider only the one-mode approximation 
in the ansatz; we did not account for the surface energy 
between the triangular phase and the vacancies; and we 
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FIG. 3: PFC simulations with different values of po for 
r = —0.9. The number of atoms increases with po. (a)-(f) 
correspond to p = 0.06, 0.08, 0.10, 0.12, 0.14 and 0.16. 

did not account for thermal fluctuations, introduced in 
the simulation to help the system equilibrate faster. 

With the modifications described above,the PFC sim- 
ulation operates very much like a molecular dynamics 
simulation, but on diffusive time scales many orders of 
magnitude faster than pure molecular dynamics We 
can control the number of atoms and the temperature in 
the system by adjusting p and the magnitude of thermal 
noise, 77, respectively. The interaction potential between 
individual PFC atoms is specified by the PFC free energy 
(specifically the gradient terms) and is controlled by the 
undercooling r. In fact, by decreasing the value of po 
such that the system is dilute enough, we can simulate 
a liquid using the PFC model! We simulated such a liq- 
uid with parameters r = —0.9, po = 0.09, a = 15 and 
(3 = 0.9. A typical result is shown in Fig. (0(b)). Fig. ([5j) 
shows the two point correlation function, g(x), extracted 
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FIG. 4: The PFC atomic density increases linearly with the 
order parameter, po, when the vacancy term is added to the 
model, r = —0.9 is used. The curve starts to saturate at 
around po — 0.15, as o pposed to the analytical prediction 
p y/(-48 - 56r)/133 0.134. 



from the simulation. It resembles the two point corre- 
lation function of a liquid — a correlation hole, a strong 
nearest neighbor correlation and a weak correlation with 
atoms one or two atomic spacings awayjlj. 
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FIG. 5: The two point correlation function of a liquid using 
the PFC model. Parameters are r = —0.9, po = 0.09, a = 15 
and p 0.9. 
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